VE_SD_SchoolConnectedness
Leigh Gayle Goetschius
Date of Last Edit: 28 January, 2021
Packages & Versions Used
# Table of packages
kable(table[-1,], format = "html", align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Package | Title | Maintainer | Version | URL | |
|---|---|---|---|---|---|
| foreign | Read Data Stored by ‘Minitab’, ‘S’, ‘SAS’, ‘SPSS’, ‘Stata’, ‘Systat’, ‘Weka’, ‘dBase’, … | R Core Team <R-core@R-project.org>; | 0.8-80 | NA | |
| tidyverse | Easily Install and Load the ‘Tidyverse’ | Hadley Wickham <hadley@rstudio.com>; | 1.3.0 | http://tidyverse.tidyverse.org, | |
| MplusAutomation | An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus | Michael Hallquist <michael.hallquist@gmail.com>; | 0.7-3 | NA | |
| interactions | Comprehensive, User-Friendly Toolkit for Probing Interactions | Jacob A. Long <long.1377@osu.edu>; | 1.1.3 | NA | |
| ggpubr | ‘ggplot2’ Based Publication Ready Plots | Alboukadel Kassambara <alboukadel.kassambara@gmail.com>; | 0.4.0 | NA | |
| psych | Procedures for Psychological, Psychometric, and Personality Research | William Revelle <revelle@northwestern.edu>; | 2.0.7 | NA | |
| table1 | Tables of Descriptive Statistics in HTML | Benjamin Rich <mail@benjaminrich.net>; | 1.2 | NA | |
| MatchIt | Nonparametric Preprocessing for Parametric Causal Inference | Kosuke Imai <kimai@Princeton.Edu>; | 3.0.2 | NA | |
| apaTables | Create American Psychological Association (APA) Style Tables | David Stanley <dstanley@uoguelph.ca>; | 2.0.5 | NA | |
| car | Companion to Applied Regression | John Fox <jfox@mcmaster.ca>; | 3.0-9 | https://r-forge.r-project.org/projects/car/, | |
| MASS | Support Functions and Datasets for Venables and Ripley’s MASS | Brian Ripley <ripley@stats.ox.ac.uk>; | 7.3-51.6 | NA | |
| data.table |
Extension of data.frame
|
Matt Dowle <mattjdowle@gmail.com>; | 1.13.0 | http://r-datatable.com, https://Rdatatable.gitlab.io/data.table, | |
| kableExtra | Construct Complex Table with ‘kable’ and Pipe Syntax | Hao Zhu <haozhu233@gmail.com>; | 1.1.0 | http://haozhu233.github.io/kableExtra/, |
Import Data
FullData = read.spss('/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/FFCW_core_ffid_transfer.sav',to.data.frame=TRUE, use.value.labels = FALSE)Preparing the Data
The composite scores for violence exposure and social deprivation are calculated exactly how they’re done in Hein et al., 2020 in SCAN.
Creating Some Initial Flags
###Creating a variable to indicate whether child lived with mother most of the time at each time point
FullData$age3cust <- ifelse(FullData$m3a2 == 1 |FullData$m3a2 == 2, "1", "0")
FullData$age5cust <- ifelse(FullData$m4a2 == 1 |FullData$m4a2 == 2, "1", "0")
FullData$age9cust <- ifelse(FullData$m5a2 == 1 |FullData$m5a2 == 2, "1", "0")###Creating a variable to indicate whether mom was currently married, living with, or romantically involved with baby’s father at age 9 (since m5c4 doesn’t seem to be recognized) ####FFCWS documentation indicates that they would assign a 1 / yes for the following responses to a4 - 1(married), 4(cohabiting or living together), 5(romantically involved but living apart), -1(refused), -2(don’t know) ####Based on conversations with Colter Mitchell about refusing itself being of interest, will instead only include1, 4, or 5
FullData$con.m5c4 = ifelse(FullData$m5a4 == 1 | FullData$m5a4 == 4 | FullData$m5a4 == 5, "1", "0")###Creating a variable to indicate whether mom was currently in any relationship at each time point
FullData$age3rel = ifelse(FullData$m3d6 == 1 |FullData$m3e2 == 1, "1", "0")
FullData$age5rel = ifelse(FullData$m4d5 == 1 |FullData$m4e2 == 1, "1", "0")
FullData$age9rel = ifelse(FullData$con.m5c4 == 1 |FullData$m5d2 == 1, "1", "0")###Creating a flag to indicate that the mom said she was not in a relationship - with BF or CP
FullData$age3.norel = ifelse(FullData$m3d6 == 2 & FullData$m3e2 == 2, "1", "0")
FullData$age5.norel = ifelse(FullData$m4d5 == 2 & FullData$m4e2 == 2, "1", "0")
FullData$age9.norel = ifelse(FullData$con.m5c4 == 2 & FullData$m5d2 == 2, "1", "0")###Creating a flag to indicate that the mom refused to answer question about relationship with BF
FullData$age3.refusebf = ifelse(FullData$m3a4 == -1, "1", "0")
FullData$age5.refusebf = ifelse(FullData$m4a4 == -1, "1", "0")
FullData$age9.refusebf = ifelse(FullData$m5a4 == -1, "1", "0")###Creating a flag to indicate that the mom refused to answer question about relationship with a CP
FullData$age3.refusecp = ifelse(FullData$m3e2 == -1, "1", "0")
FullData$age5.refusecp = ifelse(FullData$m4e2 == -1, "1", "0")
FullData$age9.refusecp = ifelse(FullData$m5d2 == -1, "1", "0")Recode for Missing
-9 (not in wave), -6 (skip), -3(missing), -2(don’t know), -1 (refuse)
is.na(FullData) <- FullData == -9
is.na(FullData) <- FullData == -8
is.na(FullData) <- FullData == -6
is.na(FullData) <- FullData == -5
is.na(FullData) <- FullData == -3
is.na(FullData) <- FullData == -2
is.na(FullData) <- FullData == -1Select Variables
Violence Exposure and Social Deprivation Composite Scores The items used in the violence exposure and social deprivation composites come from Hein et al., 2020, SCAN are are the same as those in Goetschius et al., 2020, DCN, Goetschius et al., 2020, JAMA Network Open, and Peckins et al., 2019, PNE. I haven’t repeated where they come from here, but their origins can be found in those papers. Additionally, specific papers that T.Hein selected items based on should be referenced in the code below.
School Connectedness 9 & 15 Items for this come from the Panel Study on Income Dynamics Child Development Supplement (PSID-CDS) (The Panel Study of Income Dynamics Child Development Supplement: User Guide for CDS-III, 2010)
Child Internalizing and Externalizing Symptoms Items for both of these came from the Child Behavior Checklist (Achenbach & Edelbrock, 1983). The internalizing items came from the anxious/depressed, withdrawn/depressed, and somatic symptoms subscales and the externalizing items came from the rule-breaking and aggressive behavior subscales.
Positive Adolescent Function Items from this come from the EPOCH Measure of Adolescent Wellbeing (Kern et al., 2016)
Adolescent Internalizing Symptoms Items for this come from: - Center for Epidemiologic Studies Depression Scale (CES-D) (Radloff, 1977) - Brief Symptom Inventory 18 (BSI-18) (Derogatis & Savitz, 2000)
Adolescent Externalizing Symptoms Items for this come from: - Delinquency: National Longitudinal Study of Adolescent Health (Add Health - (Harris, 2013) - Dickman’s Impulsivity scale (Dickman, 1990) - Binary substance use variables (alcohol-2+ drinks, cigarettes, marijuana, illegal drugs other than marijuana, prescription drugs w/out prescription)
Covariates Most of the covariates are pretty self-explanatory, but one requires a bit more information. To get at income I used the poverty ratio that’s constructed in the Fragile Families dataset. Because we use data from most of the waves and I don’t want to have multiple covariates representing income in my SEM for model degrees of freedom purposes, I calculate the average poverty ratio and use that as the covariate. According to the Fragile Families website, the poverty ratio is the ratio of total household income, as defined in c_1hhinc, to the official poverty thresholds, designated by the U.S. Census Bureau. The thresholds in c_1povca vary by family composition and year. At each wave, we used the poverty thresholds for the year preceding the interview. We calculated separate thresholds based on mother and father reports of household size and composition. However, calculations for married/cohabiting mothers and fathers rely on mother reports of household size and composition. A small number of missing values (don’t know, refused) were treated as 0 in household membership counts.Also included is m1city.
Primary Caregiver Relationship Used this to determine who the primary caregiver generally was for the reported data at the study waves included in the violence exposure and social deprivation composites as well as the externalizing and internalizing at age 9.
Data.ThreatDepVars = FullData %>%
dplyr::select( m1b2, cm1edu,cm2md_case_con,cm2md_case_lib,m2b17a,m2b17b,m2b17c,m2b17d,m2b17e,m2b17f,IH3j10,IH3j11,IH3j13,IH3j14,IH3j15,IH3j16,IH3j17,IH3j18,IH3j19,IH3j3,IH3j4,IH3j6,IH3j7,IH3j8,IH3j9,IH3k1a,IH3k1b,IH3k1c,IH3k1d,IH3k1e,IH3k2a,IH3k2b,IH3k2c,IH3k2d,IH3k2e,IH3l1,IH3l2,IH3l3,IH3l4,IH3l5,IH3l6,IH3l7,m3a2,m3a4,m3d6,m3d9n1,m3d9p,m3d9p1,m3d10,m3d10d,m3d1f,m3d7a,m3d7b,m3d7c,m3d7d,m3d7e,m3d7f,m3d7g,m3d7h,m3d7i,m3d7j,m3d7k,m3d7l,m3d7n,m3d7n1,m3d7p,m3d7p1,m3e2,m3e21b,m3e21f,m3e23a,m3e23b,m3e23c,m3e23d,m3e23e,m3e23f,m3e23g,m3e23h,m3e23i,m3e23j,m3e23k,m3e23l,m3e24,m3h3,m3h3a,m3h4,m3h5,m3h6,m3i6c, m3i6f,m3i23f,m3j36j,m3j36k,m3k26a, m3k26b,IH3r10,IH3s2,IH3s3,IH3s4,IH5g3,IH5g4,IH5g6,IH5g7,IH5g8,IH5g9,IH5g10,IH5g11,IH5g13,IH5g14,IH5g15,IH5g16,IH5g17,IH5g18,IH5g19,IH5h1,IH5h2,IH5h3,IH5h4,IH5h5,IH5h6,IH5h7,m4a2,m4a4,m4d1b,m4d1f,m4d5,m4d7a,m4d7b,m4d7c,m4d7d,m4d7e,m4d7f,m4d7g,m4d7h,m4d7i,m4d7j,m4d7o,m4d7p,m4d10,m4d10a,m4d10e,m4e2,m4e21b,m4e21f,m4e23a,m4e23b,m4e23c,m4e23d,m4e23e,m4e23f,m4e23g,m4e23h,m4e23i,m4e23j,m4e23o,m4e23p,m4e23q,m4e24,m4e24d,m4h3,m4h3a,m4h4,m4h5,m4h6,m4i0m1,m4i0m2,m4i0m3,m4i0m4,m4i0m5,m4i0n1,m4i0n2,m4i0n3,m4i0n4,m4i0n5,m4i23g,m4i23i,m4j22j,m4j22k,m4k25a,m4k25b,IH5r10,IH5s2,IH5s3,IH5s4,IH5s5,o5c10,o5d2,o5d3,o5d4,o5d5,n5g1f,n5g1h,hv5_dsae,hv5_dspr,hv5_dsraw,hv5_dsss,hv5_ppvtae,hv5_ppvtpr,hv5_ppvtraw,hv5_ppvtss,p5m2a,p5m2b,p5m2c,p5m2d,p5m2e,p5m3a,p5m3b,p5m3c,p5m3d,p5m3e,p5m5a,p5m5b,p5m5c,m5a2,m5a4,m5c2b,m5c2f,m5c6a,m5c6b,m5c6c,m5c6d,m5c6e,m5c6f,m5c6g,m5c6h,m5c6i,m5c6j,m5c6o,m5c6p,m5c7,m5c8a,m5c8e,m5d2,m5d19b,m5d19f,m5d20a,m5d20b,m5d20c,m5d20d,m5d20e,m5d20f,m5d20g,m5d20h,m5d20i,m5d20j,m5d20o,m5d20p,m5d21,m5d22,m5d22d,m5e3,m5e3a,m5e4,m5e5,m5e6,m5g21a,m5g21b,m5g21c,m5g21d,m5g21e,m5g21f,m5g21g,m5g21h,m5g21i,m5g21k,m5i25a,m5i25b,p5q1c,p5q1d,p5q1f,p5q1g,p5q1h,p5q1i,p5q1j,p5q1k,p5q1m,p5q1n,p5q2a,p5q2b,p5q2c,p5q2d,p5q2e,ff_id,m1city,age3cust:age9.refusecp)
Data.PAF = FullData %>%
dplyr::select("ff_id", "k6d2b", "k6d2e", "k6d2f", "k6d2g", "k6d2h", "k6d2i", "k6d2k", "k6d2l", "k6d2m", "k6d2o", "k6d2s", "k6d2u", "k6d2v", "k6d2w", "k6d2y", "k6d2aa", "k6d2ad", "k6d2ae", "k6d2af", "k6d2ah")
Data.AnxDep = FullData %>%
dplyr::select(ff_id, k6d2ag, k6d2ai, k6d2d, k6d2j, k6d2t, k6d2ac, k6d2ak, k6d2c, k6d2n, k6d2x, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66)
Data.SchoolConn = FullData %>%
dplyr::select("ff_id", "k5e1a", "k5e1b", "k5e1c", "k5e1d", "k6b1a", "k6b1b", "k6b1c", "k6b1d")
Extern_Variables = FullData %>%
dplyr::select(ff_id, k6d61a:k6d61m, k6d2a, k6d2p, k6d2r, k6d2z, k6d2ab, k6d2aj, k6d40, k6d48, k6f63, k6f68, k6f74, p6b35, p6b37:p6b39, p6b41:p6b45, p6b57, p6b59, p6b49:p6b51, p6b60:p6b64, p6b67)
CBCL = FullData %>%
dplyr::select(ff_id,p5q3m,p5q3ab,p5q3ac,p5q3ad,p5q3ae,p5q3af,p5q3ah,p5q3ar,p5q3av,p5q3ax,p5q3bq,p5q3ck,p5q3db,p5q3e,p5q3ao,p5q3bk,p5q3bo,p5q3bu,p5q3cu,p5q3cv,p5q3da,p5q3as,p5q3au,p5q3aw,p5q3az,p5q3bb1,p5q3bb2,p5q3bb3,p5q3bb4,p5q3bb5,p5q3bb6,p5q3bb7,p5q3b,p5q3x,p5q3aa,p5q3al,p5q3ap,p5q3bi,p5q3bm,p5q3br,p5q3bs,p5q3bz,p5q3ca,p5q3cj,p5q3cp,p5q3cr,p5q3ct,p5q3cx,p5q3cy,p5q3c,p5q3o,p5q3r,p5q3s,p5q3t,p5q3u,p5q3v,p5q3aj,p5q3bc,p5q3bn,p5q3cf,p5q3cg,p5q3ch,p5q3ci,p5q3cn,p5q3co,p5q3cq,p5q3cw)
Demo_Variables = FullData %>%
dplyr::select(ff_id, ck6ethrace, cm1bsex, cm1inpov, cm2povco, cm3povco, cm4povco, cm5povco, cp6povco, m1city)
PC_Info = FullData %>%
dplyr::select(ff_id, IH3int5, IH5int5, pcg5idstat)
CBCL_15 = FullData %>%
dplyr::select(ff_id, p6b36, p6b40, p6b52, p6b53, p6b54, p6b68, p6b65, p6b66, p6b35, p6b37, p6b38, p6b39, p6b41, p6b42, p6b43, p6b44, p6b45, p6b57, p6b59, p6b49, p6b50, p6b51, p6b60, p6b61, p6b62, p6b63, p6b64, p6b67)Recode Variables
School Connectedness
We need an extra NA code for the school connectedness variables (particularly the age 15 variables) – 7 on the school connectedness scale means (NA - Homeschooled)
is.na(Data.SchoolConn) <- Data.SchoolConn == 7
Data.SchoolConn$k6b1a_r = ((Data.SchoolConn$k6b1a - 5)*-1)
Data.SchoolConn$k6b1b_r = ((Data.SchoolConn$k6b1b - 5)*-1)
Data.SchoolConn$k6b1c_r = ((Data.SchoolConn$k6b1c - 5)*-1)
Data.SchoolConn$k6b1d_r = ((Data.SchoolConn$k6b1d - 5)*-1)Teen Report Internalizing
Data.AnxDep$k6d2ag_r = ((Data.AnxDep$k6d2ag-5)*-1)
Data.AnxDep$k6d2ai_r = ((Data.AnxDep$k6d2ai-5)*-1)
Data.AnxDep$k6d2d_r = ((Data.AnxDep$k6d2d-5)*-1)
Data.AnxDep$k6d2j_r = ((Data.AnxDep$k6d2j-5)*-1)
Data.AnxDep$k6d2t_r = ((Data.AnxDep$k6d2t-5)*-1)
Data.AnxDep$k6d2ac_r = ((Data.AnxDep$k6d2ac-5)*-1)
Data.AnxDep$k6d2ak_r = ((Data.AnxDep$k6d2ak-5)*-1)
Data.AnxDep$k6d2c_r = ((Data.AnxDep$k6d2c-5)*-1)
Data.AnxDep$k6d2n_r = ((Data.AnxDep$k6d2n-5)*-1)
Data.AnxDep$k6d2x_r = ((Data.AnxDep$k6d2x-5)*-1)Teen Report Externalizing
Extern_Variables = Extern_Variables %>%
mutate(k6d2a_r = (k6d2a-5)*-1,
k6d2p_r = (k6d2p-5)*-1,
k6d2r_r = (k6d2r-5)*-1,
k6d2z_r = (k6d2z-5)*-1,
k6d2ab_r = (k6d2ab-5)*-1,
k6d2aj_r = (k6d2aj-5)*-1,
k6d40_r = (k6d40-3)*-1,
k6d48_r = (k6d48-3)*-1,
k6f63_r = (k6f63-3)*-1,
k6f68_r = (k6f68-3)*-1,
k6f74_r = (k6f74-2)*-1)Positive Adolescent Function
Data.PAF$k6d2b_r = ((Data.PAF$k6d2b-5)*-1)
Data.PAF$k6d2b_r = 5-Data.PAF$k6d2b
Data.PAF$k6d2e_r = ((Data.PAF$k6d2e-5)*-1)
Data.PAF$k6d2f_r = ((Data.PAF$k6d2f-5)*-1)
Data.PAF$k6d2g_r = ((Data.PAF$k6d2g-5)*-1)
Data.PAF$k6d2h_r = ((Data.PAF$k6d2h-5)*-1)
Data.PAF$k6d2i_r = ((Data.PAF$k6d2i-5)*-1)
Data.PAF$k6d2k_r = ((Data.PAF$k6d2k-5)*-1)
Data.PAF$k6d2l_r = ((Data.PAF$k6d2l-5)*-1)
Data.PAF$k6d2m_r = ((Data.PAF$k6d2m-5)*-1)
Data.PAF$k6d2o_r = ((Data.PAF$k6d2o-5)*-1)
Data.PAF$k6d2s_r = ((Data.PAF$k6d2s-5)*-1)
Data.PAF$k6d2u_r = ((Data.PAF$k6d2u-5)*-1)
Data.PAF$k6d2v_r = ((Data.PAF$k6d2v-5)*-1)
Data.PAF$k6d2w_r = ((Data.PAF$k6d2w-5)*-1)
Data.PAF$k6d2y_r = ((Data.PAF$k6d2y-5)*-1)
Data.PAF$k6d2aa_r = ((Data.PAF$k6d2aa-5)*-1)
Data.PAF$k6d2ad_r = ((Data.PAF$k6d2ad-5)*-1)
Data.PAF$k6d2ae_r = ((Data.PAF$k6d2ae-5)*-1)
Data.PAF$k6d2af_r = ((Data.PAF$k6d2af-5)*-1)
Data.PAF$k6d2ah_r = ((Data.PAF$k6d2ah-5)*-1)Demographic Variable Prep
# Race/Ethnicity at Age 15- ck6ethrace
# 5 Multi-racial, non-hispanic
# 4 Other only, non-hispanic
# 3 Hispanic/Latino
# 2 Black/Af. American only, non-hispanic
# 1 White only, non-hispanic
Demo_Variables = Demo_Variables %>%
rowwise() %>%
mutate(povco_sum = sum(cm1inpov, cm2povco, cm3povco, cm4povco, cm5povco, cp6povco, na.rm=TRUE)) %>%
mutate(wv1_pov_inc = if_else(is.na(cm1inpov), 0,1),
wv2_pov_inc = if_else(is.na(cm2povco), 0,1),
wv3_pov_inc = if_else(is.na(cm3povco), 0,1),
wv4_pov_inc = if_else(is.na(cm4povco), 0,1),
wv5_pov_inc = if_else(is.na(cm5povco), 0,1),
wv6_pov_inc = if_else(is.na(cp6povco), 0,1)) %>%
mutate(pov_inc_sum = sum(wv1_pov_inc,wv2_pov_inc,wv3_pov_inc,wv4_pov_inc,wv5_pov_inc,wv6_pov_inc, na.rm = TRUE)) %>%
mutate(povco_avg = povco_sum/pov_inc_sum) %>%
mutate(Race_AA = if_else(ck6ethrace == 2, 1, 0),
Race_C = if_else(ck6ethrace == 1, 1, 0),
Race_L = if_else(ck6ethrace == 3, 1, 0)) %>%
mutate(cm1bsex = cm1bsex-1) %>%
dplyr::select(ff_id, povco_avg, Race_AA, Race_C, Race_L, ck6ethrace, cm1bsex, m1city)Merge to One Dataframe
All_Variables = Data.ThreatDepVars %>%
left_join(Data.AnxDep, by="ff_id")
All_Variables = All_Variables %>%
left_join(Extern_Variables, by="ff_id")
All_Variables = All_Variables %>%
left_join(Data.PAF, by = "ff_id")
All_Variables = All_Variables %>%
left_join(Data.SchoolConn, by = "ff_id")
All_Variables = All_Variables %>%
left_join(CBCL, by="ff_id")
All_Variables = All_Variables %>%
left_join(Demo_Variables, by = "ff_id")Write Variables to CSV
write.csv(All_Variables,'/Volumes/csmonk/csmonk-lab/Data/Fragile_Families_Longitudinal/Leigh_FFCWS/VE_SD_IntExtPAF_SchoolConn_012021.csv')Write Variables for MPlus
All_Variables[is.na(All_Variables)] = 99prepareMplusData(All_Variables,'All_Variables_012021.dat')Recode Missing Back to NA for Further Use
All_Variables[is.na(All_Variables)] = 99Actual Analysis Starts Here - Read in Data
Here I am reading in data written out from the creation/recoding of variables above as well as the factor scores extracted from Mplus for the simple slopes analysis. Convention for how factor scores should be extracted when probing latent variable interactions has not been clearly established. Therefore, we used two different approaches and tested whether the method used influenced the interpretation of the interaction. The first method extracted the factor scores from a measurement model containing all of the latent variables in the model (e.g. school connectedness at age 9, school connectedness at age 15, and positive adolescent function). The second approach used the factor scores extracted from individual measurement models. Due to the correlations between latent variables that are accounted for in measurement models (Kline, 2015), it is likely that the extracted factors are different depending on how they are extracted. This is why there are two dataframes, All_Variables and All_Variables2.
All_Variables = read_csv('VE_SD_IntExtPAF_SchoolConn_012021.csv')## Parsed with column specification:
## cols(
## .default = col_double()
## )## See spec(...) for full column specifications.# List of IDs that are included in models with covariates
fs_IDs = read_csv('ff_ids.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## city = col_double()
## )fs_IDs$ff_id = as.integer(fs_IDs$ff_id)
# Factor scores from SC9,SC15,PAF Measurement Mod (stratified)
fs_scores_pafmod = read_csv('factor_scores_SC15_SC9_PAF_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC15 = col_double(),
## SC15_SE = col_double(),
## SC9 = col_double(),
## SC9_SE = col_double(),
## PAF = col_double(),
## PAF_SE = col_double(),
## city = col_double()
## )fs_scores_pafmod$ff_id = as.integer(fs_scores_pafmod$ff_id)## Warning: NAs introduced by coercion# Factor scores from SC9,Ext9,Int9 Measurement Mod (stratified)
fs_scores_ext9mod = read_csv('factor_scores_SC_Int_Ext_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC15 = col_double(),
## SC15_SE = col_double(),
## SC9 = col_double(),
## SC9_SE = col_double(),
## Int15 = col_double(),
## Int15_SE = col_double(),
## Ext15 = col_double(),
## Ext15_SE = col_double(),
## Int9 = col_double(),
## Int9_SE = col_double(),
## Ext9 = col_double(),
## Ext9_SE = col_double(),
## city = col_double()
## )fs_scores_ext9mod$ff_id = as.integer(fs_scores_ext9mod$ff_id)## Warning: NAs introduced by coercion# Individual factor scores -- needs to be updated to stratified
SC9 = read_csv('factor_scores_SC9_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## SC9 = col_double(),
## SC9_SE = col_double(),
## city = col_double()
## )SC9$ff_id = as.integer(SC9$ff_id)## Warning: NAs introduced by coercionSC15 = read_csv('factor_scores_SC15_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## SC15 = col_double(),
## SC15_SE = col_double(),
## city = col_double()
## )SC15$ff_id = as.integer(SC15$ff_id)
Int9 = read_csv('factor_scores_Int9_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## INT9 = col_double(),
## INT9_SE = col_double(),
## city = col_double()
## )Int9$ff_id = as.integer(Int9$ff_id)## Warning: NAs introduced by coercionExt9 = read_csv('factor_scores_Ext9_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_character(),
## EXT9 = col_double(),
## EXT9_SE = col_double(),
## city = col_double()
## )Ext9$ff_id = as.integer(Ext9$ff_id)## Warning: NAs introduced by coercionInt15 = read_csv('factor_scores_Int15_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## INT15 = col_double(),
## INT15_SE = col_double(),
## city = col_double()
## )Int15$ff_id = as.integer(Int15$ff_id)
Ext15 = read_csv('factor_scores_Ext15_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## EXT15 = col_double(),
## EXT15_SE = col_double(),
## city = col_double()
## )Ext15$ff_id = as.integer(Ext15$ff_id)
PAF = read_csv('factor_scores_PAF_012221.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## PAF = col_double(),
## PAF_SE = col_double(),
## city = col_double()
## )PAF$ff_id = as.integer(PAF$ff_id)
# Merge things
PAF_mod_data = fs_IDs %>%
left_join(fs_scores_pafmod) %>%
left_join(All_Variables) ## Joining, by = c("ff_id", "city")## Joining, by = "ff_id"Ext9_mod_data = fs_IDs %>%
left_join(fs_scores_ext9mod) %>%
left_join(All_Variables) ## Joining, by = c("ff_id", "city")
## Joining, by = "ff_id"Indiv_FS_All_Variables = fs_IDs %>%
left_join(All_Variables) %>%
left_join(SC9) %>%
left_join(SC15) %>%
left_join(Int9) %>%
left_join(Ext9) %>%
left_join(Int15) %>%
left_join(Ext15) %>%
left_join(PAF)## Joining, by = "ff_id"## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")
## Joining, by = c("ff_id", "city")This reverses any 99s that are marked in the MPlus data – this should already be done, but just in case.
is.na(PAF_mod_data) = PAF_mod_data == 99
is.na(Ext9_mod_data) = Ext9_mod_data == 99
is.na(Indiv_FS_All_Variables) = Indiv_FS_All_Variables == 99Look at distributions of our variables
PAF = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=PAF)) + xlab("Positive Adolescent Function") + geom_density() + theme_classic()
Dep = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=DepComp)) + xlab("Social Deprivation") + geom_density() + theme_classic()
Threat = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=ThreatComp)) + xlab("Violence Exposure") + geom_density() + theme_classic()
SC9 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=SC9)) + xlab("School Connectedness Age 9") + geom_density() + theme_classic()
SC15 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=SC15)) + xlab("School Connectedness Age 15") + geom_density() + theme_classic()
Int9 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=INT9)) + xlab("Internalizing Psychopathology Age 9") + geom_density() + theme_classic()
Ext9 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=EXT9)) + xlab("Externalizing Psychopathology Age 9") + geom_density() + theme_classic()
Int15 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=INT15)) + xlab("Internalizing Psychopathology Age 15") + geom_density() + theme_classic()
Ext15 = Indiv_FS_All_Variables %>%
ggplot(mapping = aes(x=EXT15)) + xlab("Externalizing Psychopathology Age 15") + geom_density() + theme_classic()
figure <- ggarrange(PAF, Threat, Dep, SC9,SC15, Int9, Ext9, Int15,Ext15,
ncol = 2, nrow = 5)## Warning: Removed 2 rows containing non-finite values (stat_density).## Warning: Removed 353 rows containing non-finite values (stat_density).## Warning: Removed 53 rows containing non-finite values (stat_density).## Warning: Removed 366 rows containing non-finite values (stat_density).
## Warning: Removed 366 rows containing non-finite values (stat_density).## Warning: Removed 2 rows containing non-finite values (stat_density).## Warning: Removed 1 rows containing non-finite values (stat_density).figure
Skewness and Kurtosis
psych::describe(Indiv_FS_All_Variables$ThreatComp) # -- this is potentially problematic, it seems like there are some serious kurtosis issues ## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0.01 0.53 -0.09 -0.05 0.45 -1.13 7.1 8.23 1.95 11.84 0.01psych::describe(Indiv_FS_All_Variables$DepComp)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0 0.53 -0.08 -0.04 0.47 -1.47 4.02 5.49 1.34 4.3 0.01psych::describe(Indiv_FS_All_Variables$PAF)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 3244 -0.02 0.89 -0.01 -0.01 0.93 -3.04 1.56 4.6 -0.15 -0.39
## se
## X1 0.02psych::describe(Indiv_FS_All_Variables$INT9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 2880 0.08 0.82 0.03 0.02 0.86 -1.01 4.92 5.92 0.73 1.38 0.02psych::describe(Indiv_FS_All_Variables$EXT9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 2880 0.05 0.88 0.02 0 0.95 -1.21 4.28 5.48 0.44 0.06 0.02psych::describe(Indiv_FS_All_Variables$INT15)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3244 0.03 0.87 0.02 0 0.95 -1.37 3.18 4.55 0.24 -0.46 0.02psych::describe(Indiv_FS_All_Variables$EXT15)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3245 0.04 0.91 0.07 0.06 0.9 -2.06 2.8 4.86 -0.13 -0.22 0.02psych::describe(Indiv_FS_All_Variables$SC9)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 2893 -0.05 0.77 0.01 0.02 0.9 -2.28 0.86 3.14 -0.59 -0.2 0.01psych::describe(Indiv_FS_All_Variables$SC15)## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 3193 -0.04 0.78 -0.02 0.03 0.76 -2.65 0.91 3.57 -0.56 -0.22
## se
## X1 0.01Table Creation
The code below creates: - Demographic Table - Correlation Table w/descriptives for our main variables.
# Import a categorical race variable because I had accidentally cut it out earlier.
Race_Cat = read_csv('Race_Categorical.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## Race_Cat = col_character()
## )Race_Cat$Race_Cat = as.factor(Race_Cat$Race_Cat)
MatBirthVars = read_csv('Maternal_Birth_Variables.csv')## Parsed with column specification:
## cols(
## ff_id = col_double(),
## m1b2 = col_double(),
## cm1edu = col_double()
## )MatBirthVars$m1b2 = as.factor(MatBirthVars$m1b2)
MatBirthVars$cm1edu = as.factor(MatBirthVars$cm1edu)
Indiv_FS_All_Variables = Indiv_FS_All_Variables %>%
left_join(Race_Cat, by='ff_id') %>%
left_join(MatBirthVars, by='ff_id')
Indiv_FS_All_Variables = Indiv_FS_All_Variables %>%
mutate(MatMar = case_when(
m1b2 == 1 ~ 'Married',
m1b2 == 2 ~ 'Not Married',
TRUE ~ 'NA'
),
MatEdu = case_when(
cm1edu == 1 ~ 'Less than high school',
cm1edu == 2 ~ 'High school or equivalent',
cm1edu == 3 ~ 'Some college or technical school',
cm1edu == 4 ~ 'College or graduate school',
TRUE ~ 'NA'
),
sex = case_when(
cm1bsex == 0 ~ 'Male',
cm1bsex == 1 ~ 'Female',
TRUE ~ 'NA'
))table1::label(Indiv_FS_All_Variables$sex) = "Sex"
table1::label(Indiv_FS_All_Variables$ck6ethrace) = "Race-Ethnicity"
table1::label(Indiv_FS_All_Variables$povco_avg) = "Average Poverty Ratio"
table1::label(Indiv_FS_All_Variables$MatMar) = "Maternal Marital Status at Child's Birth"
table1::label(Indiv_FS_All_Variables$MatEdu) = "Maternal Education at Child's Birth"
table1::label(Indiv_FS_All_Variables$ThreatComp) = "Violence Exposure"
table1::label(Indiv_FS_All_Variables$DepComp) = "Social Deprivation"
table1::label(Indiv_FS_All_Variables$SC9) = "School Connectedness Age 9"
table1::label(Indiv_FS_All_Variables$SC15) = "School Connectedness Age 15"
table1::label(Indiv_FS_All_Variables$PAF) = "Positive Adolescent Function"
table1::label(Indiv_FS_All_Variables$INT15) = "Internalizing Symptoms Age 15"
table1::label(Indiv_FS_All_Variables$EXT15) = "Externalizing Symptoms Age 15"
table1::label(Indiv_FS_All_Variables$INT9) = "Internalizing Symptoms Age 9"
table1::label(Indiv_FS_All_Variables$EXT9) = "Externalizing Symptoms Age 9"
table1::table1(~sex + Race_Cat + povco_avg + MatMar + MatEdu, data = Indiv_FS_All_Variables)| Overall (N=3246) |
|
|---|---|
| Sex | |
| Female | 1585 (48.8%) |
| Male | 1661 (51.2%) |
| Race_Cat | |
| African American | 1592 (49.0%) |
| Caucasian | 587 (18.1%) |
| Latinx | 808 (24.9%) |
| Other | 259 (8.0%) |
| Average Poverty Ratio | |
| Mean (SD) | 2.11 (2.10) |
| Median [Min, Max] | 1.46 [0.120, 21.2] |
| Maternal Marital Status at Child's Birth | |
| Married | 785 (24.2%) |
| NA | 18 (0.6%) |
| Not Married | 2443 (75.3%) |
| Maternal Education at Child's Birth | |
| College or graduate school | 365 (11.2%) |
| High school or equivalent | 1030 (31.7%) |
| Less than high school | 1025 (31.6%) |
| NA | 5 (0.2%) |
| Some college or technical school | 821 (25.3%) |
table1::table1(~ThreatComp + DepComp + SC9 + SC15 + INT9 + EXT9 + INT15 + EXT15 + PAF, data = Indiv_FS_All_Variables)| Overall (N=3246) |
|
|---|---|
| Violence Exposure | |
| Mean (SD) | 0.00869 (0.535) |
| Median [Min, Max] | -0.0922 [-1.13, 7.10] |
| Social Deprivation | |
| Mean (SD) | 0.00323 (0.533) |
| Median [Min, Max] | -0.0780 [-1.47, 4.02] |
| School Connectedness Age 9 | |
| Mean (SD) | -0.0524 (0.769) |
| Median [Min, Max] | 0.0100 [-2.28, 0.865] |
| Missing | 353 (10.9%) |
| School Connectedness Age 15 | |
| Mean (SD) | -0.0440 (0.784) |
| Median [Min, Max] | -0.0190 [-2.66, 0.912] |
| Missing | 53 (1.6%) |
| Internalizing Symptoms Age 9 | |
| Mean (SD) | 0.0752 (0.821) |
| Median [Min, Max] | 0.0310 [-1.01, 4.92] |
| Missing | 366 (11.3%) |
| Externalizing Symptoms Age 9 | |
| Mean (SD) | 0.0467 (0.878) |
| Median [Min, Max] | 0.0230 [-1.21, 4.28] |
| Missing | 366 (11.3%) |
| Internalizing Symptoms Age 15 | |
| Mean (SD) | 0.0257 (0.873) |
| Median [Min, Max] | 0.0185 [-1.37, 3.18] |
| Missing | 2 (0.1%) |
| Externalizing Symptoms Age 15 | |
| Mean (SD) | 0.0418 (0.908) |
| Median [Min, Max] | 0.0710 [-2.06, 2.80] |
| Missing | 1 (0.0%) |
| Positive Adolescent Function | |
| Mean (SD) | -0.0182 (0.894) |
| Median [Min, Max] | -0.0140 [-3.04, 1.56] |
| Missing | 2 (0.1%) |
CorTable = Indiv_FS_All_Variables %>%
dplyr::select(ThreatComp, DepComp, SC9, SC15, INT9, EXT9, INT15, EXT15, PAF)
apa.cor.table(CorTable, filename = '/Users/leighgoetschius/Box Sync/Shift_Persist_FF/Manuscript/Figures/cor_table_012321.doc', table.number = NA,
show.conf.interval = TRUE, landscape = TRUE)Quick check of VE-SD associations
cor.test(Indiv_FS_All_Variables$ThreatComp, Indiv_FS_All_Variables$DepComp)##
## Pearson's product-moment correlation
##
## data: Indiv_FS_All_Variables$ThreatComp and Indiv_FS_All_Variables$DepComp
## t = 24.925, df = 3244, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3716339 0.4293926
## sample estimates:
## cor
## 0.4009115fit_cor = lm(PAF ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor2 = lm(EXT9 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor3 = lm(INT9 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor4 = lm(EXT15 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
fit_cor5 = lm(INT15 ~ ThreatComp + DepComp, data = Indiv_FS_All_Variables)
paf_vif = car::vif(fit_cor)
ext9_vif = car::vif(fit_cor2)
int9_vif = car::vif(fit_cor3)
ext15_vif = car::vif(fit_cor4)
int15_vif = car::vif(fit_cor5)
mean(paf_vif[1],ext9_vif[1],int9_vif[1],ext15_vif[1],int15_vif[1])## [1] 1.190807Plot of the SC15 Effect
PAF_form = lm(Indiv_FS_All_Variables$PAF~Indiv_FS_All_Variables$SC15)
Int_form = lm(Indiv_FS_All_Variables$INT15~Indiv_FS_All_Variables$SC15)
Ext_form = lm(Indiv_FS_All_Variables$EXT15~Indiv_FS_All_Variables$SC15)
PAF_SC15 = Indiv_FS_All_Variables %>%
ggplot(aes(x=SC15, y=PAF)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="skyblue3") +
xlab("School Connectedness Age 15") +
ylab("Positive Adolescent Function Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = PAF_form) +
theme_classic() +
theme(text = element_text(family='serif'))
Int_SC15 = Indiv_FS_All_Variables %>%
ggplot(aes(x=SC15, y=INT15)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="red") +
xlab("School Connectedness Age 15") +
ylab("Internalizing Symptoms Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int_form) +
theme_classic() +
theme(text = element_text(family='serif'))
Ext_SC15 = Indiv_FS_All_Variables %>%
ggplot(aes(x=SC15, y=EXT15)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour = "red")+
xlab("School Connectedness Age 15") +
ylab("Externalizing Symptoms Age 15") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext_form) +
theme_classic() +
theme(text = element_text(family='serif'))
SC15_figure <- ggarrange(Int_SC15, Ext_SC15, PAF_SC15,
ncol = 2, nrow = 2)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 55 rows containing non-finite values (stat_smooth).## Warning: Removed 55 rows containing non-finite values (stat_regline_equation).## Warning: Removed 55 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 54 rows containing non-finite values (stat_smooth).## Warning: Removed 54 rows containing non-finite values (stat_regline_equation).## Warning: Removed 54 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 55 rows containing non-finite values (stat_smooth).## Warning: Removed 55 rows containing non-finite values (stat_regline_equation).## Warning: Removed 55 rows containing missing values (geom_point).SC15_figure
Plot SC9 Contemporaneous Effect
Int9_form = lm(Indiv_FS_All_Variables$INT9~Indiv_FS_All_Variables$SC9)
Ext9_form = lm(Indiv_FS_All_Variables$EXT9~Indiv_FS_All_Variables$SC9)
Int9_SC9 = Indiv_FS_All_Variables %>%
ggplot(aes(x=SC9, y=INT9)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour="red") +
xlab("School Connectedness Age 9") +
ylab("Internalizing Symptoms Age 9") +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Int9_form) +
theme_classic() +
theme(text = element_text(size=20, family='serif'))
Ext9_SC9 = Indiv_FS_All_Variables %>%
ggplot(aes(x=SC9, y=EXT9)) +
geom_point(position="jitter") +
geom_smooth(method="lm", colour = "red")+
xlab("School Connectedness Age 9") +
ylab("Externalizing Symptoms Age 9") +
theme_classic() +
theme(text = element_text(size=20, family='serif')) +
stat_regline_equation(aes(label = paste(..adj.rr.label..)),formula = Ext9_form)
SC9_figure <- ggarrange(Int9_SC9, Ext9_SC9,
ncol = 2, nrow = 1)## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 399 rows containing non-finite values (stat_smooth).## Warning: Removed 399 rows containing non-finite values (stat_regline_equation).## Warning: Removed 399 rows containing missing values (geom_point).## `geom_smooth()` using formula 'y ~ x'## Warning: Removed 399 rows containing non-finite values (stat_smooth).## Warning: Removed 399 rows containing non-finite values (stat_regline_equation).## Warning: Removed 399 rows containing missing values (geom_point).SC9_figure
Probing SC9xSD –> PAF w.stratification
This uses factor scores extracted from a measurement model with SC9, SC15, and PAF that is stratified by m1city (city at birth)
SC9 as moderator
# Center/Standardize the Variables
PAF_mod_data = PAF_mod_data %>%
mutate(PAF_scaled = scale(PAF),
DepComp_scaled = scale(DepComp),
SC9_scaled = scale(SC9))
# Set your model to a variable
sc9xpaf_mod = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data)
sc9PAF_ros = data.frame(xmin = -2.772,
xmax = 2.766,
ymin = -Inf,
ymax = Inf)
# Interpret the interaction
interact_plot(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "Social Deprivation", y.label = "Positive Function Age 15", legend.main = "School Connectedness", data = PAF_mod_data) + geom_rect(data=sc9PAF_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.3, inherit.aes=FALSE) + geom_vline(xintercept = 0, colour="black", linetype = "longdash") + theme_classic() + theme(text=element_text(size = 16, family="serif"))
johnson_neyman(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-15.914, -1.270], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.262, 1.456]
johnson_neyman(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [2.766, 32.588], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
sim_slopes(sc9xpaf_mod, pred = DepComp_scaled, modx = SC9_scaled, data = PAF_mod_data, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-15.914, -1.270], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.262, 1.456]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -1.000 (- 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.066 0.024 -2.729 0.006
##
## Slope of DepComp_scaled when SC9_scaled = -0.000 (Mean):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.109 0.017 -6.381 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.000 (+ 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.151 0.025 -5.949 0.000SD as moderator
SDPAF_ros = data.frame(xmin = -1.270,
xmax = 1.456,
ymin = -Inf,
ymax = Inf)
interact_plot(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "School Connectedness (Age 9)", y.label = "Positive Function Age 15", legend.main = "Social Deprivation", data = PAF_mod_data) + geom_rect(data=SDPAF_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.2, inherit.aes=FALSE) + geom_vline(xintercept = 0, colour="black", linetype = "longdash") + theme_classic() + theme(text=element_text(size = 20, family="serif"))
sim_slopes(sc9xpaf_mod, pred = SC9_scaled, modx = DepComp_scaled, data = PAF_mod_data, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [2.766, 32.588], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of SC9_scaled when DepComp_scaled = -1.000 (- 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.263 0.025 10.636 0.000
##
## Slope of SC9_scaled when DepComp_scaled = -0.000 (Mean):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.221 0.017 12.981 0.000
##
## Slope of SC9_scaled when DepComp_scaled = 1.000 (+ 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.179 0.025 7.160 0.000Getting coefficients, variances, and covariances for and running the Preacher interaction tool.
vcov(sc9xpaf_mod)## (Intercept) DepComp_scaled SC9_scaled
## (Intercept) 2.896838e-04 8.515363e-07 1.622477e-07
## DepComp_scaled 8.515363e-07 2.901365e-04 1.892200e-05
## SC9_scaled 1.622477e-07 1.892200e-05 2.896314e-04
## DepComp_scaled:SC9_scaled 2.129388e-05 1.310645e-05 2.497242e-06
## DepComp_scaled:SC9_scaled
## (Intercept) 2.129388e-05
## DepComp_scaled 1.310645e-05
## SC9_scaled 2.497242e-06
## DepComp_scaled:SC9_scaled 3.277457e-04summary(sc9xpaf_mod)##
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3807 -0.6616 0.0105 0.6985 2.2747
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002748 0.017020 -0.161 0.8717
## DepComp_scaled -0.108682 0.017033 -6.381 2.02e-10 ***
## SC9_scaled 0.220912 0.017019 12.981 < 2e-16 ***
## DepComp_scaled:SC9_scaled -0.042298 0.018104 -2.336 0.0195 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9674 on 3242 degrees of freedom
## Multiple R-squared: 0.06504, Adjusted R-squared: 0.06418
## F-statistic: 75.18 on 3 and 3242 DF, p-value: < 2.2e-16psych::describe(PAF_mod_data$DepComp_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0 1 -0.15 -0.09 0.89 -2.77 7.54 10.31 1.34 4.3 0.02psych::describe(PAF_mod_data$SC9_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0 1 0.07 0.08 1 -3.26 1.46 4.72 -0.63 0.11 0.02Preacher Output - SC9 as Moderator
Your Input
X1 = -2.77 X2 = 7.54 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = -0.002748 X Slope = -0.108682 Z Slope = 0.220912 XZ Slope = -0.042298 df = 3240 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.00028968 var(b1) 0.00029014 var(b2) 0.00028963 var(b3) 0.00032775 cov(b2,b0) 1.6e-7 cov(b3,b1) 0.00001311
Region of Significance
Z at lower bound of region = -15.914 Z at upper bound of region = -1.2702 (simple slopes are significant outside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = -0.2237(0.0241), t=-9.2951, p=0 simple slope = -0.0664(0.0243), t=-2.7291, p=0.0064 At Z = cv2… simple intercept = -0.0027(0.017), t=-0.1615, p=0.8717 simple slope = -0.1087(0.017), t=-6.3805, p=0 At Z = cv3… simple intercept = 0.2182(0.0241), t=9.0616, p=0 simple slope = -0.151(0.0254), t=-5.949, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = -3.5183(0.2714), t=-12.9657, p=0 simple slope = 0.5644(0.2879), t=1.9607, p=0.05 Upper Bound…
simple intercept = -0.2834(0.0275), t=-10.3015, p=0 simple slope = -0.055(0.028), t=-1.9606, p=0.05
Points to Plot
Line for cv1: From {X=-2.77, Y=-0.0398} to {X=7.54, Y=-0.7242} Line for cv2: From {X=-2.77, Y=0.2983} to {X=7.54, Y=-0.8222} Line for cv3: From {X=-2.77, Y=0.6364} to {X=7.54, Y=-0.9202}
Preacher Simple Slopes Plot
xx <- c(-2.77,7.54) # <-- change to alter plot dims
yy <- c(-1.2315,0.6364) # <-- change to alter plot dims
leg <- c(-2.77,-0.9981) # <-- change to alter legend location
x <- c(-2.77,7.54) # <-- x-coords for lines
y1 <- c(-0.0398,-0.7242)
y2 <- c(0.2983,-0.8222)
y3 <- c(0.6364,-0.9202)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='X',ylab='Y',main='MLR 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))
Testing if SDxSC9 Interaction Improves Model Fit
I calculated this based on this MPlus page: https://www.statmodel.com/chidiff.shtml
https://stats.idre.ucla.edu/mplus/faq/how-can-i-compute-a-chi-square-test-for-nested-models-with-the-mlr-or-mlm-estimators/#:~:text=Typically%20a%20chi%2Dsquare%20difference,freedom%20between%20the%20two%20models.
It references this paper: Satorra, A., & Bentler, P. M. (2010). Ensuring Positiveness of the Scaled Difference Chi-square Test Statistic. Psychometrika, 75(2), 243–248. https://doi.org/10.1007/s11336-009-9135-y
To Compute the Satorra-Bentler Scaled Chi-Square Difference Test:
To calculate the scaling factor:
cd=(p0 * c0 - p1*c1)/(p0 - p1) Where p=parameters and c=correction factor
Then, the test statistic is as follows TRd = -2*(L0 - L1)/cd Where L=loglikelihood
X0 values are the ones with fewer freely estimated parameters X1 values are for the moderation model
This is compared to the chi-square distribution based on the degrees of freedom difference between the two models, in this case, it is 2. Therefore, our model with the interactions fits better than the main effects models.
Compared to ME Only Model
p0=97
p1=99
c0=1.0271
c1=1.0268
cd = (p0*c0-p1*c1)/(p0-p1)
cd## [1] 1.01225L0=-50172.965
L1=-50169.544
TRd = -2*(L0-L1)/cd
TRd## [1] 6.7592df_paf=p1-p0
qchisq(.95, df_paf)## [1] 5.991465pchisq(TRd, df=df_paf, lower.tail=FALSE)## [1] 0.03406108Compared to Model Where SDxSC9 Int. Set to 0
p0=98
p1=99
c0=1.0270
c1=1.0268
cd = (p0*c0-p1*c1)/(p0-p1)
cd## [1] 1.0072L0=-50172.381
L1=-50169.544
TRd = -2*(L0-L1)/cd
TRd## [1] 5.633439df_paf=p1-p0
qchisq(.95, df_paf)## [1] 3.841459pchisq(TRd, df=df_paf, lower.tail=FALSE)## [1] 0.01762103Who does this interaction not apply to?
According to the Johnson-Neyman intervals, when SC9 <-1.429 stdev and when social deprivation was > 3.067 stdev, they were outside of the range of the interaction. Let’s see who that applies to.
PAF_mod_data %>%
filter(SC9_scaled < -1.270) %>%
summarize(n(),
mean(DepComp_scaled, na.rm = TRUE),
min(DepComp_scaled),
max(DepComp_scaled),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled),
max(PAF_scaled))## # A tibble: 1 x 7
## `n()` `mean(DepComp_s… `min(DepComp_sc… `max(DepComp_sc… `mean(PAF_scale…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 364 0.0926 -1.71 4.68 -0.342
## # … with 2 more variables: `min(PAF_scaled)` <dbl>, `max(PAF_scaled)` <dbl>PAF_mod_data %>%
filter(DepComp_scaled > 2.766) %>%
summarize(n(),
mean(SC9_scaled, na.rm = TRUE),
min(SC9_scaled, na.rm = TRUE),
max(SC9_scaled, na.rm=TRUE),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled, na.rm = TRUE),
max(PAF_scaled, na.rm = TRUE))## # A tibble: 1 x 7
## `n()` `mean(SC9_scale… `min(SC9_scaled… `max(SC9_scaled… `mean(PAF_scale…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 43 -0.140 -3.15 1.29 -0.322
## # … with 2 more variables: `min(PAF_scaled, na.rm = TRUE)` <dbl>,
## # `max(PAF_scaled, na.rm = TRUE)` <dbl>PAF_mod_data %>%
filter(DepComp_scaled > 2.776 & SC9_scaled < -1.270) %>%
summarize(n(),
mean(PAF_scaled, na.rm = TRUE),
min(PAF_scaled, na.rm = TRUE),
max(PAF_scaled, na.rm = TRUE))## # A tibble: 1 x 4
## `n()` `mean(PAF_scaled, na.rm… `min(PAF_scaled, na.rm… `max(PAF_scaled, na.rm…
## <int> <dbl> <dbl> <dbl>
## 1 3 -0.954 -1.51 0.0494Checking model assumptions
(Approximate since model is done in Mplus)
# Checking some model assumptions
mean(sc9xpaf_mod$residuals)## [1] -2.142609e-17t.test(sc9xpaf_mod$residuals)##
## One Sample t-test
##
## data: sc9xpaf_mod$residuals
## t = -1.2625e-15, df = 3245, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.03327605 0.03327605
## sample estimates:
## mean of x
## -2.142609e-17plot(sc9xpaf_mod)
sred = studres(sc9xpaf_mod)
hist(sred, freq=FALSE,
main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit)
Robustness Check of the Interaction
We tested if the method for extracting the factor scores made a difference in the interpretation of the interaction. The analyses below used factor scores extracted from CFAs containing ONLY one latent variable at a time. I’m using this as the robustness check rather than the preferred method because I think the overall latent variable moderation model accounts for correlations between the latent variables and these are picked up when extracting factor scores from CFAs with all of the latent variables included in the moderation model.
Indiv_FS_All_Variables = Indiv_FS_All_Variables %>%
mutate(PAF_scaled = scale(PAF),
DepComp_scaled = scale(DepComp),
SC9_scaled = scale(SC9))
sc9xpaf_mod_indiv = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)
summary(sc9xpaf_mod_indiv)##
## Call:
## lm(formula = PAF_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4193 -0.6713 -0.0051 0.7105 2.1402
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0006067 0.0182564 -0.033 0.9735
## DepComp_scaled -0.1107125 0.0187281 -5.912 3.79e-09 ***
## SC9_scaled 0.1220655 0.0182720 6.680 2.85e-11 ***
## DepComp_scaled:SC9_scaled -0.0436768 0.0194845 -2.242 0.0251 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9798 on 2887 degrees of freedom
## (355 observations deleted due to missingness)
## Multiple R-squared: 0.02972, Adjusted R-squared: 0.02871
## F-statistic: 29.47 on 3 and 2887 DF, p-value: < 2.2e-16johnson_neyman(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-20.12, -1.21], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.89, 1.19]
johnson_neyman(sc9xpaf_mod_indiv, pred = SC9_scaled, modx = DepComp_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is OUTSIDE the interval [1.35, 22.44], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-1.91, 7.54]
sim_slopes(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, data = Indiv_FS_All_Variables, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is OUTSIDE the interval [-20.124, -1.210], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.894, 1.193]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -0.999 (- 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.067 0.026 -2.536 0.011
##
## Slope of DepComp_scaled when SC9_scaled = 0.001 (Mean):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.111 0.019 -5.913 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.000 (+ 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.154 0.028 -5.602 0.000interact_plot(sc9xpaf_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = Indiv_FS_All_Variables) + theme_classic() + theme(text=element_text(size = 16, family="serif"))
Probing SC9xSD –> Ext9 w.stratification
Ext9_mod_data = Ext9_mod_data %>%
mutate(DepComp_scaled = scale(DepComp),
SC9_scaled = scale(SC9),
Ext9_scaled = scale(Ext9))
sc9xext9_mod = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data)sc9ext9_ros = data.frame(xmin = -2.9,
xmax = 3.221,
ymin = -Inf,
ymax = Inf)
interact_plot(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, interval = FALSE, rug = TRUE, rug.sides = 'bl', modx.labels = c("Low**","Mean**","High**"), x.label = "Social Deprivation", y.label = "Externalizing Symptoms Age 9", legend.main = "School Connectedness", data = Ext9_mod_data) + theme_classic() + theme(text=element_text(size = 20, family="serif")) + geom_rect(data=sc9ext9_ros, aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax),fill="lightgray", linetype=0, alpha=.3, inherit.aes=FALSE) + geom_vline(xintercept = 0, colour="black", linetype = "longdash") 
interact_plot(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, interval = TRUE, x.label = "School Connectedness", y.label = "Externalizing Sx Age 9", legend.main = "Social Deprivation", data = Ext9_mod_data) + theme_classic() 
johnson_neyman(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is INSIDE the interval [-3.325, 106.882], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.300, 1.595]
johnson_neyman(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is INSIDE the interval [-105.418, 3.221], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
sim_slopes(sc9xext9_mod, pred = DepComp_scaled, modx = SC9_scaled, data = Ext9_mod_data, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is INSIDE the interval [-3.325, 106.882], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-3.300, 1.595]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -1.000 (- 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.194 0.024 8.191 0.000
##
## Slope of DepComp_scaled when SC9_scaled = -0.000 (Mean):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.227 0.017 13.617 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.000 (+ 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.259 0.025 10.424 0.000sim_slopes(sc9xext9_mod, pred = SC9_scaled, modx = DepComp_scaled, data = Ext9_mod_data, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is INSIDE the interval [-105.418, 3.221], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-2.772, 7.540]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of SC9_scaled when DepComp_scaled = -1.000 (- 1 SD):
##
## Est. S.E. t val. p
## -------- ------- --------- -------
## -0.254 0.024 -10.529 0.000
##
## Slope of SC9_scaled when DepComp_scaled = -0.000 (Mean):
##
## Est. S.E. t val. p
## -------- ------- --------- -------
## -0.222 0.017 -13.331 0.000
##
## Slope of SC9_scaled when DepComp_scaled = 1.000 (+ 1 SD):
##
## Est. S.E. t val. p
## -------- ------- -------- -------
## -0.189 0.024 -7.746 0.000Getting coefficients, variances, and covariances for and running the Preacher interaction tool.
vcov(sc9xext9_mod)## (Intercept) DepComp_scaled SC9_scaled
## (Intercept) 2.767544e-04 1.144093e-06 2.565935e-07
## DepComp_scaled 1.144093e-06 2.772599e-04 2.227941e-05
## SC9_scaled 2.565935e-07 2.227941e-05 2.766394e-04
## DepComp_scaled:SC9_scaled 2.504654e-05 1.430275e-05 3.207774e-06
## DepComp_scaled:SC9_scaled
## (Intercept) 2.504654e-05
## DepComp_scaled 1.430275e-05
## SC9_scaled 3.207774e-06
## DepComp_scaled:SC9_scaled 3.131165e-04summary(sc9xext9_mod)##
## Call:
## lm(formula = Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5730 -0.6417 -0.0311 0.5608 5.8561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002607 0.016636 0.157 0.8755
## DepComp_scaled 0.226741 0.016651 13.617 <2e-16 ***
## SC9_scaled -0.221720 0.016632 -13.331 <2e-16 ***
## DepComp_scaled:SC9_scaled 0.032589 0.017695 1.842 0.0656 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9444 on 3242 degrees of freedom
## Multiple R-squared: 0.109, Adjusted R-squared: 0.1082
## F-statistic: 132.2 on 3 and 3242 DF, p-value: < 2.2e-16psych::describe(Ext9_mod_data$SC9_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0 1 0.09 0.08 1.04 -3.3 1.59 4.89 -0.64 0.14 0.02psych::describe(Ext9_mod_data$DepComp_scaled)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 3246 0 1 -0.15 -0.09 0.89 -2.77 7.54 10.31 1.34 4.3 0.02Preacher Output - SC9 as Moderator
Your Input
X1 = -2.77 X2 = 7.54 cv1 = -1 cv2 = 0 cv3 = 1 Intercept = 0.002607 X Slope = 0.226741 Z Slope = -0.22172 XZ Slope = 0.032589 df = 3240 alpha = 0.05
Asymptotic (Co)variances
var(b0) 0.00027675 var(b1) 0.00027726 var(b2) 0.00027664 var(b3) 0.00031312 cov(b2,b0) 2.6e-7 cov(b3,b1) 0.0000143
Region of Significance
Z at lower bound of region = -3.3254 Z at upper bound of region = 106.8632 (simple slopes are significant inside this region.)
Simple Intercepts and Slopes at Conditional Values of Z
At Z = cv1… simple intercept = 0.2243(0.0235), t=9.5404, p=0 simple slope = 0.1942(0.0237), t=8.1915, p=0 At Z = cv2… simple intercept = 0.0026(0.0166), t=0.1567, p=0.8755 simple slope = 0.2267(0.0167), t=13.6172, p=0 At Z = cv3… simple intercept = -0.2191(0.0235), t=-9.31, p=0 simple slope = 0.2593(0.0249), t=10.4235, p=0
Simple Intercepts and Slopes at Region Boundaries
Lower Bound…
simple intercept = 0.7399(0.0577), t=12.814, p=0 simple slope = 0.1184(0.0604), t=1.9607, p=0.05 Upper Bound…
simple intercept = -23.6911(1.7775), t=-13.3284, p=0 simple slope = 3.7093(1.8918), t=1.9607, p=0.05
Points to Plot
Line for cv1: From {X=-2.77, Y=-0.3135} to {X=7.54, Y=1.6882} Line for cv2: From {X=-2.77, Y=-0.6255} to {X=7.54, Y=1.7122} Line for cv3: From {X=-2.77, Y=-0.9375} to {X=7.54, Y=1.7362}
Preacher Simple Slopes Plot
xx <- c(-2.77,7.54) # <-- change to alter plot dims
yy <- c(-1.4722,1.7362) # <-- change to alter plot dims
leg <- c(-2.77,-1.0711) # <-- change to alter legend location
x <- c(-2.77,7.54) # <-- x-coords for lines
y1 <- c(-0.3135,1.6882)
y2 <- c(-0.6255,1.7122)
y3 <- c(-0.9375,1.7362)
plot(xx,yy,type='n',font=2,font.lab=2,xlab='X',ylab='Y',main='MLR 2-Way Interaction Plot')
lines(x,y1,lwd=3,lty=1,col=1)
lines(x,y2,lwd=3,lty=5,col=2)
lines(x,y3,lwd=3,lty=6,col=3)
points(x,y1,col=1,pch=16)
points(x,y2,col=1,pch=16)
points(x,y3,col=1,pch=16)
legend(leg[1],leg[2],legend=c('CVz1(1)','CVz1(2)','CVz1(3)'),lwd=c(3,3,3),lty=c(1,5,6),col=c(1,2,3))
Testing if the SDxSC9 Interaction Improves Fit
Compared to ME Only Model
p0=187
p1=191
c0=1.0505
c1=1.0345
cd = (p0*c0-p1*c1)/(p0-p1)
cd## [1] 0.2865L0=-70459.229
L1=-70454.076
TRd = -2*(L0-L1)/cd
TRd## [1] 35.97208df_paf=p1-p0
qchisq(.95, df_paf)## [1] 9.487729pchisq(TRd, df=df_paf, lower.tail=FALSE)## [1] 2.932224e-07Compared to Model Where SDxSC9 Int. Set to 0
p0=190
p1=191
c0=1.0386
c1=1.0345
cd = (p0*c0-p1*c1)/(p0-p1)
cd## [1] 0.2555L0=-70456.106
L1=-70454.076
TRd = -2*(L0-L1)/cd
TRd## [1] 15.89041df_paf=p1-p0
qchisq(.95, df_paf)## [1] 3.841459pchisq(TRd, df=df_paf, lower.tail=FALSE)## [1] 6.711791e-05Who does this interaction not apply to?
Ext9_mod_data %>%
filter(DepComp_scaled > 2.493) %>%
summarize(n(),
mean(SC9_scaled, na.rm = TRUE),
min(SC9_scaled, na.rm = TRUE),
max(SC9_scaled, na.rm=TRUE),
mean(Ext9_scaled, na.rm = TRUE),
min(Ext9_scaled, na.rm = TRUE),
max(Ext9_scaled, na.rm = TRUE))## # A tibble: 1 x 7
## `n()` `mean(SC9_scale… `min(SC9_scaled… `max(SC9_scaled… `mean(Ext9_scal…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 58 -0.191 -3.14 1.28 0.506
## # … with 2 more variables: `min(Ext9_scaled, na.rm = TRUE)` <dbl>,
## # `max(Ext9_scaled, na.rm = TRUE)` <dbl>Ext9_mod_data %>%
filter(SC9_scaled < -2.749) %>%
summarize(n(),
mean(DepComp_scaled, na.rm = TRUE),
min(DepComp_scaled, na.rm = TRUE),
max(DepComp_scaled, na.rm=TRUE),
mean(Ext9_scaled, na.rm = TRUE),
min(Ext9_scaled, na.rm = TRUE),
max(Ext9_scaled, na.rm = TRUE))## # A tibble: 1 x 7
## `n()` `mean(DepComp_s… `min(DepComp_sc… `max(DepComp_sc… `mean(Ext9_scal…
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 40 0.0132 -1.49 2.81 0.657
## # … with 2 more variables: `min(Ext9_scaled, na.rm = TRUE)` <dbl>,
## # `max(Ext9_scaled, na.rm = TRUE)` <dbl>Ext9_mod_data %>%
filter(DepComp_scaled > 2.493 & SC9_scaled < -2.749) %>%
summarize(n(),
mean(Ext9_scaled, na.rm = TRUE),
min(Ext9_scaled, na.rm = TRUE),
max(Ext9_scaled, na.rm = TRUE))## # A tibble: 1 x 4
## `n()` `mean(Ext9_scaled, na.r… `min(Ext9_scaled, na.r… `max(Ext9_scaled, na.r…
## <int> <dbl> <dbl> <dbl>
## 1 1 0.801 0.801 0.801Checking model assumptions
(Approximate since model is done in Mplus)
# Checking some model assumptions
mean(sc9xext9_mod$residuals)## [1] -5.269587e-18t.test(sc9xext9_mod$residuals)##
## One Sample t-test
##
## data: sc9xext9_mod$residuals
## t = -3.1806e-16, df = 3245, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.0324847 0.0324847
## sample estimates:
## mean of x
## -5.269587e-18plot(sc9xext9_mod)
sred = studres(sc9xext9_mod)
hist(sred, freq=FALSE,
main="Distribution of Studentized Residuals")
xfit<-seq(min(sred),max(sred),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit)
Robustness Check – Alternative Way of Extracting Factor Scores
Indiv_FS_All_Variables = Indiv_FS_All_Variables %>%
mutate(Ext9_scaled = scale(EXT9))
sc9xext9_mod_indiv = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Indiv_FS_All_Variables)
johnson_neyman(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is INSIDE the interval [-3.27, 99.13], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.89, 1.19]
johnson_neyman(sc9xext9_mod_indiv, pred = SC9_scaled, modx = DepComp_scaled)## JOHNSON-NEYMAN INTERVAL
##
## When DepComp_scaled is INSIDE the interval [-56.59, 1.77], the slope of
## SC9_scaled is p < .05.
##
## Note: The range of observed values of DepComp_scaled is [-1.91, 7.54]
sim_slopes(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, data = Indiv_FS_All_Variables, digits = 3)## JOHNSON-NEYMAN INTERVAL
##
## When SC9_scaled is INSIDE the interval [-3.273, 99.125], the slope of
## DepComp_scaled is p < .05.
##
## Note: The range of observed values of SC9_scaled is [-2.894, 1.193]
##
## SIMPLE SLOPES ANALYSIS
##
## Slope of DepComp_scaled when SC9_scaled = -0.999 (- 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.207 0.026 7.928 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 0.002 (Mean):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.242 0.018 13.114 0.000
##
## Slope of DepComp_scaled when SC9_scaled = 1.002 (+ 1 SD):
##
## Est. S.E. t val. p
## ------- ------- -------- -------
## 0.278 0.027 10.219 0.000interact_plot(sc9xext9_mod_indiv, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Externalizing Symptoms (Age 9)", legend.main = "School Connectedness", data = Indiv_FS_All_Variables) + theme_classic() + theme(text=element_text(size = 16, family="serif"))
Hand plotting interaction
PAF_mod_data = PAF_mod_data %>%
mutate(SC_bin = case_when(
SC9_scaled <= -1 ~ 'Low',
SC9_scaled <= 1 ~ 'Med',
SC9_scaled > 1 ~ 'High',
TRUE ~ 'NA'
))
Ext9_mod_data = Ext9_mod_data %>%
mutate(SC_bin = case_when(
SC9_scaled <= -1 ~ 'Low',
SC9_scaled <= 1 ~ 'Med',
SC9_scaled > 1 ~ 'High',
TRUE ~ 'NA'
))
PAF_mod_data_bin = PAF_mod_data %>%
filter(SC_bin != 'NA')
Ext9_mod_data_bin = Ext9_mod_data %>%
filter(SC_bin != 'NA')PAF_mod_data_bin %>%
ggplot(aes(x=DepComp_scaled, y=PAF_scaled, group=SC_bin, color=SC_bin)) +
stat_smooth(se = TRUE) +
xlab("Social Deprivation (Centered)") + ylab("Positive Function (Centered)")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
PAF_mod_data_bin %>%
ggplot(aes(x=DepComp_scaled, y=PAF_scaled, group=SC_bin, color=SC_bin)) +
stat_smooth(method = "lm", se = TRUE) +
xlab("Social Deprivation (Centered)") + ylab("Positive Function (Centered)")## `geom_smooth()` using formula 'y ~ x'
interact_mod_PAF = lm(PAF_scaled ~ DepComp_scaled * SC9_scaled, data = PAF_mod_data_bin)
interact_plot(interact_mod_PAF, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation", y.label = "Positive Adolescent Function", legend.main = "School Connectedness", data = PAF_mod_data_bin) + theme_classic() + theme(text=element_text(size = 16, family="serif"))
Ext9_mod_data_bin %>%
ggplot(aes(x=DepComp_scaled, y=Ext9_scaled, group=SC_bin, color=SC_bin)) +
stat_smooth(se = TRUE) +
xlab("Social Deprivation (Centered)") + ylab("Externalizing Age 9 (Centered)")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Ext9_mod_data_bin %>%
ggplot(aes(x=DepComp_scaled, y=Ext9_scaled, group=SC_bin, color=SC_bin)) +
stat_smooth(method="lm",se = TRUE) +
xlab("Social Deprivation (Centered)") + ylab("Externalizing Age 9 (Centered)")## `geom_smooth()` using formula 'y ~ x'
interact_mod_Ext9 = lm(Ext9_scaled ~ DepComp_scaled * SC9_scaled, data = Ext9_mod_data_bin)
interact_plot(interact_mod_Ext9, pred = DepComp_scaled, modx = SC9_scaled, interval = TRUE, modx.labels = c("Low","Mean","High"), x.label = "Social Deprivation (Centered)", y.label = "Externalizing Age 9 (Centered)", legend.main = "School Connectedness", data = Ext9_mod_data_bin) + theme_classic() + theme(text=element_text(size = 16, family="serif"))
Who is the Primary Caregiver?
Age 3: IH3int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather
Age5: IH5int5 01 Bio Mother 02 Bio Father 03 Mat Grandmother 04 Pat Grandmother 05 Other Relative 06 Other Non-Relative 07 Mat Grandmother 08 Mat Grandfather
Age 9: pcg5idstat 61 Bio Mother
63 Other
62 Bio Father
PC_Info = read_csv('primary_caregiver_relationships.csv')## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_double(),
## ff_id = col_double(),
## IH3int5 = col_double(),
## IH5int5 = col_double(),
## pcg5idstat = col_double()
## )PC_Info %>%
group_by(IH3int5) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 7 x 2
## IH3int5 `n()`
## <dbl> <int>
## 1 1 3191
## 2 2 15
## 3 3 10
## 4 4 3
## 5 5 4
## 6 6 5
## 7 NA 1670PC_Info %>%
group_by(IH5int5) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 7 x 2
## IH5int5 `n()`
## <dbl> <int>
## 1 1 2914
## 2 2 27
## 3 3 19
## 4 4 7
## 5 5 10
## 6 6 1
## 7 NA 1920PC_Info %>%
group_by(pcg5idstat) %>%
summarize(n())## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 4 x 2
## pcg5idstat `n()`
## <dbl> <int>
## 1 61 3828
## 2 62 180
## 3 63 194
## 4 NA 696Hypothetical Resilience Graphs
Adapted from Luthar et al., 2000
Risk = c(2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4,2,4)
Competence = c(2,1,3,2,4,3,
2,1,3,1.5,5,2.5,
4,2,4.5,3,5,5,
2,1,3,3,4,5)
Moderator = c("Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High","Low","Low","Medium","Medium","High","High")
Type = c("promotive","promotive","promotive","promotive","promotive","promotive","prot-reac","prot-reac","prot-reac","prot-reac","prot-reac","prot-reac","prot-stab","prot-stab","prot-stab","prot-stab","prot-stab","prot-stab", "prot-enh", "prot-enh", "prot-enh", "prot-enh", "prot-enh", "prot-enh")
hypothetical_plots = data.frame(Risk, Competence, Moderator, Type)
hypothetical_plots$Moderator <- factor(hypothetical_plots$Moderator, levels = c("Low","Medium","High"))P = hypothetical_plots %>%
filter(Type == "promotive") %>%
ggplot(aes(Risk, Competence, group=Moderator)) +
geom_line(aes(linetype=Moderator), color="black") +
scale_x_continuous(limits = c(1.9,4.1)) +
scale_y_continuous(limits = c(0.9,4.1)) +
ggtitle("Promotive") +
theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face="bold"), legend.position="none", text=element_text(family="serif"))
PE = hypothetical_plots %>%
filter(Type == "prot-enh") %>%
ggplot(aes(Risk, Competence, group=Moderator)) +
geom_line(aes(linetype=Moderator), color="black") +
scale_x_continuous(limits = c(1.9,4.1)) +
scale_y_continuous(limits = c(0.9,5.1)) +
ggtitle("Protective-Enhancing") +
theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), legend.position="none", text=element_text(family="serif"))
PR = hypothetical_plots %>%
filter(Type == "prot-reac") %>%
ggplot(aes(Risk, Competence, group=Moderator)) +
geom_line(aes(linetype=Moderator), color="black") +
scale_x_continuous(limits = c(1.9,4.1)) +
scale_y_continuous(limits = c(0.9,5.1)) +
ggtitle("Protective-Reactive") +
theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), legend.position="none", text=element_text(family="serif"))
PS = hypothetical_plots %>%
filter(Type == "prot-stab") %>%
ggplot(aes(Risk, Competence, group=Moderator)) +
geom_line(aes(linetype=Moderator), color="black") +
scale_x_continuous(limits = c(1.9,4.1)) +
scale_y_continuous(limits = c(0.9,5.1)) +
ggtitle("Protective-Stabilizing") +
theme(axis.ticks = element_blank(), axis.text = element_blank(), panel.background = element_blank(), panel.grid = element_blank(), axis.line = element_line(), plot.title = element_text(hjust = 0.5, face = "bold"), text=element_text(family="serif"), legend.text=element_text(size = 12), legend.position = "bottom", legend.title = element_text(size=12, face="bold"), legend.key = element_rect(fill = "white")) +
guides(linetype = guide_legend(title.position="top", title.hjust = 0.5), override.aes = list())Hyp_figure <- ggarrange(P, PE, PR, PS, common.legend = TRUE,legend.grob = get_legend(PS), legend="bottom",
ncol = 2, nrow = 2)
Hyp_figure 